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In many biological processes heterogeneity within cell populations is an important 
issue. In this work we consider populations where the behavior of every single 
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Q cell can be described by a system of ordinary differential equations. Heterogeneity 

among individual cells is accounted for by differences in parameter values and 
initial conditions. Hereby, parameter values and initial conditions are subject to 

^ a distribution function which is part of the model specification. Based on the 

single cell model and the considered parameter distribution, a partial differential 
equation model describing the distribution of cells in the state and in the output 
space is derived. 

For the estimation of the parameter distribution within the model, we consider 
experimental data as obtained from flow cytometric analysis. From these noise- 
corrupted data a density-based statistical data model is derived. Using this data 
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model the parameter distribution within the cell population is computed using 
convex optimization techniques. 

To evaluate the proposed method, a model for the caspase activation cascade 
is considered. It is shown that for known noise properties the unknown parameter 
distributions in this model are well estimated by the proposed method. 

Keywords: parameter estimation, cell population, kernel-density estimation, flow 
cytometry, convex optimization 



1 Introduction 



Most of the modeling performed in the area of systems biology aims at achieving 
a quantitative description of intracellular pathways. Hence, most available models 
describe a "typical cell" on the basis of experimental data. Unfortunately, experi- 
mental data are in general obtained using cell population experiments, e.g. western 
blotting. If the considered population is highly heterogeneous, meaning that there 
is a large cell-cell variability, fitting a single cell model to cell population data can 
lead to biologically meaningless results. To understand the dynamical behavior of 
heterogeneous cell populations it is crucial to develop integrated cell population 
models. 

Modeling on the population scale has already been addressed by |Mantzaris| 
(2007) and Munsky et al. (2009). These authors demonstrated that populations 



can show a bimodal response if stochasticity in biochemical reactions is considered. 
But besides stochasticity in biochemical reactions there are other reasons which 
can also lead to heterogeneity in populations. Examples are unequal partition- 
ing of cellular material at cell division (Mantzaris, 2007), genetic and epigenetic 



differences (Avery, 2006). 



For the purpose of this paper, we describe heterogeneity in populations by dif- 
ferences in parameter values of the model describing the single cell dynamics. The 
network structure is assumed to be identical in all cells, as this usually represents 
the physical interactions among molecules, which should be independent of the 
cell's state. This parametric approach is well suited for genetic and epigenetic 
differences. The distribution of parameter values within the cell population of in- 
terest is described by a multivariate probability density function, which is part of 
the model specification. 

In the following the problem of estimating the parameter distribution func- 
tion is studied. Therefore, we consider high-throughput experimental methods 
such as flow cytometry, which can be used to measure concentration distributions 
within cell populations by suitable fluorescent labeled antibodies. Classical flow 
cytometry devices can measure several thousand cells per second. 
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To estimate the parameter distributions, in a first step, an appropriate pop- 
ulation model has to be found. In the literature mathematical models of cell 



populations are either described as cell ensembles (Waldherr et al.l 2009 Munsky 



et al. 2009), or as a non-linear partial differential equation (PDE) for the distri- 



bution of the state variables ( 


Mantzaris , 


2007; 


Luzyanina et al. 2009 


Tsuchiya 



et al. , 1966). In case of ensemble models, a differential equation is assigned to 
each cell, making an in depth theoretical analysis difficult. PDE models, which 
describe the time evolution of the distributions of the state variables based on the 
single cell models, are easy to handle from a theoretical point of view but hard 
to simulate for a large state dimension of the single cell model. Therefore, only 
low dimensional PDE models of populations have been studied in literature so far 



(Mantzaris, 2007 Luzyanina et al., 2007, 2009). 



In this paper a PDE model for the state distribution within a heterogeneous 
cell population is derived. Given the solution of this PDE the probability density 
of measuring a certain output can be determined. As for the estimation only the 
measured outputs are required, a numerical method for computing the output dis- 



tribution is outlined. This methods employs a particle-based approach (Rawlings 



and Bakshi, 2006) and classical density estimation (Silverman 1986). 

Based on these efficient computation scheme for the population response an 
estimation method is developed. A statistical model of the measured output dis- 
tribution is derived from the single cell measurement obtained at every measure- 
ment instance. Therefore, again kernel density estimators are used as they have 



better asymptotic properties than commonly used naive estimators (Luzyanina 



et al. 2009). Given a model and the output distribution estimated from the mea- 
surement, a /2-norm minimization is performed over the set of possible parameter 
distributions. By employing the model properties and a parameterization of the 
parameter distribution this optimization problem is convex and can be solved ef- 
ficiently. 

The paper is structured as follows. In Section 2, the problem of estimating the 
parameter distribution is introduced. In Section 3, we present the statistical model 
for the measured data and the simulation model for state and output distribution. 
Section 4 gives a short overview of the used identification procedure before in 
Section 5 the proposed method is applied to a caspase activation model with 
artificial data. 

Notation: Consider the m-dimensional hypersurface S C W 1 . The integral / 
of a function f(x), with / : W 1 — > R, over x G S is written as 



f(x)dS. 



Furthermore, the i.th unit vector is denoted by ej. 
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2 Problem statement 



For the purpose of this work, a model of a biochemical reaction network in a 
population of M cells is given by the collection of differential equations 

y ® = h{x {i \p {i) ), ie{l,...,M} 

with state variables x®(t) £ W 1 , measured variables y^\t) £ M m , and parameters 
pW £ y§q_ The index i specifies the individual cells within the population. The 
parameters can be kinetic constants, e.g. reaction rates or binding affinities. 
The cell-cell interaction of the considered pathway is assumed to be negligible, as 
it is the case in many in vitro lab experiments. 

In the following heterogeneity within the cell population is introduced, modeled 
by differential parameter values and initial conditions among individual cells. The 
distribution of parameters anc l initial conditions Xq is given by a probability 
density function $ : M. n+q — > M + with f Rn+q §(xo,p)dxodp = 1. For ease of nota- 
tion, we write £o = (^o\p T ) T - The probability density function $ is part of the 
model specification and the parameters and initial conditions of cell i are subject 
to the probability distribution 

Pr(eg <&,•", < Zn +q ) = f 1 --- Hi)dii ■ ■ ■ d£ n+q . (3) 

As outlined in Section [TJ for the study of cell populations high-throughput 
cell population measurements are available. Using these experimental techniques 
protein concentrations within thousands of cells can be measured at every mea- 
surement instance, tk, k — 1, . . . , N. This yields the measurement data 

V k = {(t k ,^(t k ))}. &k , k = l,...,N (4) 

where ip^ is the measured output of the cell i and X k is the index set of the 
cells measured at time t k . Note that the cells cannot be tracked over time, and 
are removed from the population in order to obtain the measurements. Thus, 
no single-cell time series data are available. On the other hand, the samples are 
independent and equally distributed and card(Xfc) is assumed to be large, such 
that an approximation of the output distribution is possible. 

Like most measurement devices, also high-throughput fluorescence measure- 
ments are subject to noise. For the rest of the paper, noise consisting of a relative 
and an absolute part is considered, 

^ i \t k ) = dmg( V 1 )y^(t k ) + r ] 2 , (5) 
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in which ip^ is the measured output and if e IR m is a vector of log-normally 
distributed random variables with probability density functions 

yielding the joint probability density 

m 

©w)=n e K# ( 7 ) 

0=1 

Log-normally distributed random variables are chosen here, since they are a good 
model for the commonly seen noise distributions of the considered measurement 
device and conserve the positivity of all variables. For notational simplicity the 
measurement errors of the different concentrations are assumed to be uncorrelated. 
This constraint can be removed easily. 

Given this setup the problem we are concerned with is: 

Problem 1 Given the measurement data T>k, k = 1,...,N, the cell population 
model and the noise model (j7]) ; determine the parameter distribution $(£). 

Unfortunately, estimation of using a cell population model with a finite 
number of cells and discrete sampled data is fairly difficult as no single cell tra- 
jectories are available. A far more natural approach would be to use a density de- 
scription, as the available measurement data can be interpreted as samples drawn 
from the probability density function of the output. This interpretation is also 
quite appealing from a point of modeling as the number of cells considered in a 
standard lab experiment is of the order of 10 9 and hence nevertheless too large 
to be simulated on an individual basis. In the next chapter a PDE model for the 
probability density of the output and a density model for the measurement data 
is derived. 

3 Density-based modeling of heterogeneous cell 
populations 

As outlined in the previous section, a continuous statistical model for the mea- 
surement data, as well as for the evolution of the state and output density would 
be preferable. These two aspects are addressed in the following. 
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3.1 Density model of measurement data 



The data collected by the considered measurement devices T>k are samples drawn 
from the distribution of the measured output, as mentioned in Section [2] Let 
tk) be the distribution of the measured outputs ^ (*fc) & t time tk- As ^(ip, tk) 
is considered to be a probability density, classical density estimation methods can 
be employed for estimating ^(ip, tk) from the given samples T>k- 

In this work, the problem of determining ^>(ip,tk) from T>k is approached 
using kernel density estimators. Kernel density estimators are non-parametric 



approaches to estimate probability distributions from sampled data (Silverman 



1986). They are widely used and can be thought of as placing probability "bumps" 



at each observation, as depicted in Figure [TJ These "bumps" are the kernel func- 
tion K, with L m K{ip)dip = 1. Note that here only the equations for the one 
dimensional case are given. The extension towards higher dimensions is straight- 



forward and can be found in Silverman (1986). In this work, a Gaussian kernel 
given by 

1 f 1 

2 



K(iJ)-il>®,h) 



2nh 



cxp 



h 



with standard deviation h is used. In this context, h is also called smoothing 



parameter in the literature (Silverman, 1986). 



Given the kernel K an estimator of the probability density for a given set of 
samples T>k is 



(9) 



;ex fc 



where is the cardinality of X^. The selection of the smoothing parameter h 
is crucial and depends strongly on M}~. In this work h is chosen according to the 
least-squares cross-validation method (Stone, 1984). As M\. is considered to be of 



order 10 4 , it can be assumed that the the estimated output distribution ins close 
to the actual output distribution. 



3.2 PDE model of density evolution 

As outlined previously, a continuous model for the output density is desirable 
for the purpose of parameter identification. Therefore, a PDE model for the cell 
population is derived in the next step. 

At first the single cell model is transformed in an extended state space model 

_(*«.«»>). em-{S) 
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Figure 1: Gaussian kernel density estimate ( — ) of ^(ip, t) for the measured outputs 
(o) and the associated Gaussian kernels ( — ). 

in which the parameters are appended to the state vector, £W = [^ 1,l \^ 2 '^] T G 
M. n+q with = a?W and ^ 2 ' 1 ^ = pW. This system can also be written as 

to which we refer as the extended state space representation. 



Based on ( 11 ), the PDE model for the population is derived. The state variable 
of this PDE is the state distribution function H : M n+9 X R —»• R+ :(£,£) S(f , t), 
which is defined on the extended state space. Based on the distribution function 
E, the probability of picking at random a cell from the population with states 
€ X at time t is given by 

Pr(£»(f)e*)= / H(f,*)df. (12) 

To determine the PDE for E, an infinitesimal volume = A^ i x . . . x X^ n+q of the 
extended state space is considered, with X^ = ^ + A£j]. For the 2-dimensional 
case this is depicted in Figure [2j 

For this infinitesimal volume the flux and storage balance is, 

r r JL r t+At / \ 

E({,t + M)d£- E(i,t)di=J2 / {Eitt,T)-~;^ T ))dT. (13) 

The left hand side of the equation represents the storage term and the right hand 
side the fluxes across the boundaries. The fluxes Ef and Ej are given by the 
surface integral of the product of boundary distribution and entering velocity, 
determined by the single cell dynamics, 



Et(£,t)= [ Fi(i)S(lt)dS 

JS+(£.i) 



(14) 



(£,*)= / F t (OE^t)dS, 
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Figure 2: Infinitesimal volume element X of the extended state space, with fluxes 
across the boundaries. 



in which Ss(£, i) = {£l& = 6 A li E Xtfj ^ i}. 

Next, ( 13 ) and ( flip are used to derive the PDE for the time evolution of £). 
Therefore, at first the storage term is expanded using its Taylor series, yielding 

j E&t+At)d£- I s(e l ode = (s(e,«+At)-H(e,o)II A ^ +0 ^ +l )- 

JXj: J X% j = 1 

(15) 

Here it is assumed that 0(A£j) = 0(A£) Wj G {1, . . . ,n + q}. In a second step 
the flux difference AHj(£, r) = S^~(£, r) — r) is rewritten, 



AS,«,t) = - 



<9(F,H 



A& + 0(A£))dS 



A? 



(16) 



i=i 



The first line follows from the definition of ASj(^,t) and the Taylor series expan- 
sion of Fi(£ + e;A£j)H(£ + ejA£j,t). To obtain the second line the integration is 



carried out. The final reformulation is the expansion of the time integral in (13) 
resulting in 



t+At 



N 



+ 0(A£ N )0(At 2 ). 

(17) 



Substituting (15) and (17) in the flux balance (13) and dividing by Atf]^ A£ 
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then yields, 



5(g,t + At)-5(g,t) + Q(Afl 
At 



N 

E 

i=l 



dm) 



0(At). 



Given this the PDE governing the evolution of E(£, t) is obtained by taking the 
limits A£j — > and At — > 0, leading to 



TV 



dt 



8=1 



96 



(19) 



for sufficiently smooth S(£, t). This final equation is somehow what we expected, 
a transport equation with position dependent transport direction and velocity, 



according to the single cell dynamics. The initial condition of (22) is the initial 
distribution on the extended state space, 



<(£,o) = $(0, V£e 



!>n+q 



(20) 



From the state distribution S, the output distribution T is computed as the 
integral of the state distribution along H(£) = y, 



T(i/,t) 



KU)dS, 



(21) 



where S T (y) = {£|#(0 = y}. 

The resulting partial differential equation system is 



ot 
r(y,t) 



N 

E 

i=l 



«9(f,: 



96 



-(e,t), s(e,o) = $(o 



(22) 



s(e,t)^, 



where S : R n+q x R ->■ R + and T : M m x R -> R + . This PDE is of first order, 
quasilinear and known as Liouville's equation. The solution always exists for 



sufficiently smooth F(-) (Evans, 1998) 



As the measurements are noise corrupted, the distribution of measured outputs 
ty(ip,t) is different from the actual output distribution T(y,t). It is defined by 



*(iM)=/ r^eV)© 2 ^ 2 )^, 



(23) 



where S$(ip) = {[y T , (?? 1 ) T , {rf) T ] T \d\a.g{r] l )y + if = ip}. 
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3.3 Numerical solution of PDE 



In order to study the time evolution of the output distribution T(y,t) and the 



measured output distribution equation (22) has to be solved for given <E>. 

As E(£, t) is defined on the (n + g)-dimensional space, standard grid based solvers 



are not able to solve (22 ) for n+q > 3. Theoretically, the methods of characteristics 



can be used (Evans, 1998) but for the high dimensional system we are going to 



study, also this method is difficult to apply. Instead, a stochastic method is used, 



which is known from particle filtering (Rawlings and Bakshi, 2006). 



This stochastic integration method is based on a particle description of the 
model, which is in our case equivalent to the cell ensemble model (|2]). To compute 
ty(ijj } t k ), at first a set of samples {(x^,^)}^...^, is drawn from $(£)) where S is 
the number of samples. For this set of samples the single cell model ^ is simulated, 
resulting in a set of simulated outputs {(|/ ( ' J - ) (t))}i=i ) ...,5. y^\t) is then corrupted 
by noise according to ^ resulting in {C0^(£))}i=i,...,s- Given this a numerical 
approximation of \I/(?/>,t) can be determined using the kernel density estimator 
described in Section |3.1| This numerical stochastic approximation the output of 



( 22 ) can be shown to converge as S — > oo. Hence, the measured output distribution 



\I/(?/>,tfc) can be axproximated also for high dimensional nonlinear systems. 



4 Estimation of parameter distributions 

As mentioned in Section [2] the problem studied in this work is the estimation of 
the parameter distribution $ from the data T>k- This problem is approached in 
the following by minimizing the ^-norm of the model-data mismatch, 



N 
k=l 



(24) 



in which t, $) is the distribution of the measured output ip obtained by simu- 
lation with the parameter distribution $(£). According to the cost J, the optimal 
parameter distribution $*(£) is than given by 

<f>* = argmin,| J($) 

subject to J R n+ q = 1 (25) 

HO > 0V£ G 

where the last two constraints enforce that is a probability distribution. 

Remark 1 In the whole section the measured outputs ^(tp,t) are compared with 
the noise corrupted simulated output ty(ip,t,Q). This is possible as we assume a 
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Figure 3: Schematic of head functions A*(p). 



large number of measured cells per measurement instant and therefore have good 
statistics on the measurement error. 



Unfortunately, the optimization problem (25) is infinite dimensional. There- 
fore, a parametrization of $, 



(26) 



i=l 



with a weighting vector (p e W 1 ^ is introduced. In this work the ansatz functions 
A* for <f> are chosen to be classical head functions, as depicted in Figure [3J This 
yields the simplified, finite-dimensional optimization problem, 



(f* = argmin^ J($ v ) 
subject to c T ip = 1 
<p>0, 



(27) 



in which Cj = J Rn+q A % (£)d£. The two constraints are again needed to ensure that 
is a probability density. 



In order to solve (27) using computational techniques the quasi-linearity of (22 ) 



is employed. As the superposition principle holds, the output ^(Vs t, $<p) can be 
written as the weighted sum 



(2f 



i=i 



where ^f(i/),t,A' 1 ) is the output distribution obtained for simulation with a param- 
eter distribution according to A l (£). This allows the reformulation of the objective 
function to 

N 



J $ 



fc=i 



tfftMfc) -^^^(^,t fe ,A 4 ) 



(29) 
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Employing this the optimization problem (27) can finally be written as 



Lp* = argmin^ J2k=i ( A kV ~ h) T W {A k (p - b k ) 

subject to c T f = 1 (30) 

where the integral || • \ \\ is approximated, e.g. using the trapezoidal rule. The 
column vector b k contains hereby the values ^(ip, t k ) at the grid points of the dis- 
cretization. Equivalently, the ith column of A k contains the values of ^(ip, t k , A*) 
at the grid points. The matrix W is a constant weighting matrix, determined by 
the chosen approximation of || • |||. 

Note that problem (30) is convex. Hence, even in the case of high dimensional 
if, convergence to the optimal parameter distribution within the considered class 
of distributions can be guaranteed. 



5 Application to the caspase cascade 

Programmed cell death, also called apoptosis, is an important physiological process 
to remove infected, malfunctioning, or no longer needed cells from a multicellular 
organism. Pathways to induce apoptosis converge at the caspase activation cascade 



( Hengartner 


2000 


)■ 


Eissing et al. 


( 


2004) 



Here, we consider the caspase activation in response to an 
external death receptor stimulus, e.g. the tumor necrosis factor (TNF). As seen 
from experimental cytotoxicity assays, the cellular response to a TNF stimulus is 
highly heterogeneous, with some cells dying and others surviving. To understand 
the process at the physiological level it is thus crucial to consider the cellular 
heterogeneity, using for example cell population modeling. 
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The reactions for the single cell model are given by 



C3 + C8* 




C3* + C8* 


C3* + C8 




C3* + C8* 


C3* + IAP 


k?, 

k-3 


C3* ~ IAP 


Lo + lAr 


k4 
—f 




C8* 


k^ 





C3 


kg 

-4 





C3* ~ IAP 


k^ 





IAP 


ks 





C8 


kg 
k-g 





C3 


fcio 





C8* + BAR 


fell 
fe-11 


C8* ~ BAR 


BAR 


fel2 
fe-12 





C8* ~ BAR 


fcl| 





TNFR + C8 


kij 


TNFR + C8* 



For nominal parameter values, we refer to the original publication (Eissing et al 



2004). In comparison to the original model, we added reaction vu for the initiator 
caspase 8 (C8) activation by the TNF receptor complexes (TNFR). The reaction 
rate for this activation is given by t>i 4 = £44 [TNFR] [C8], with the parameter value 
ki4 = 10~ 6 (moleculesmin) _1 . A sketch of the single cell model is given in Figure[4j 
Heterogeneity is modeled by a log-normally distributed production rate of the 
inhibitor of apoptosis IAP, kg, and a log-normally distributed amount of TNF- 
receptor complexes on the cell membrane, TNFR. These two quantities were chosen 
as it is known from experiments that there is a high cell-cell variability. Especially 
the concentration of IAPs contained in a cells is highly variable, and a variation in 



IAP production is known to affect cell death considerably (Eissing et al. , |2006[ ). In 



the following the possibility of estimating the distributions of Q(k$) and $([TNFR]) 
from the distributions of [C3*], ^(C3*), is studied. The statistical model of the 
distribution, \I/(C3*) is shown in Figure [5] This statistical model has been derived 
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Figure 4: Schematic of the caspase activation cascade. 




Figure 5: Artificial noisy measurement data for amount of active caspase 3, [C3*]. 

using artificial measurement data of 10 4 cells at the measurement instances tk, 
k = 1, . . . , 6. This is a realistic number for standard cytofluorometric experiments. 
The noise properties are assumed to be known and have been set to fi\ = 0, 
oi = 0.1, /i 2 = log(10 3 ), and <t 2 = 0.3, corresponding to an average measurement 
error of more than 20 percent. 

Based on these data, the approach presented in Section [4] is used to obtain 
an estimate for the parameter distribution. For this purpose the considered pa- 
rameter set is divided using a 12 x 12 grid, with logarithmically distributed grid 
points. The grid points are used as edge and center points of the ansatz functions 
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Figure 6: Real ( — ) vs. estimated (-o) parameter distribution, with grid points 
(o). 

A l (k 8 , [TNFR]) for $(A; 8 , [TNFR]). The obtained estimation result is depicted in 
Figure |6j 

It is obvious that the estimated parameter distribution approximates the real 
parameter distribution very well, especially considering the finite number of de- 
grees of freedom. Hence, even though there is an average measurement error of 20 
% on the single cell measurement, due to good statistics at the population level, 
the actual parameter distributions can be estimated accurately. Furthermore, this 
study shows that in principle, measuring one concentration can give enough infor- 
mation to estimate several parameter distributions, if the output distribution is 
sensitive with respect to these parameters. 

6 Summary and Conclusion 

Heterogeneity in cell populations is an important issue for research in systems 
biology. However, so far only few models describing heterogeneous populations of 
cells with more than one state variable have been developed. In this paper a partial 
differential equation model describing the time evolution of the state distribution 
is derived. We focused hereby in particular on the distribution of the measured 
outputs. 

In the second part of the paper, the model of the noise corrupted measured out- 
puts and its particular properties are used to estimate the parameter distributions 
underlying the heterogeneity. Therefore, a density-based statistical model of the 
sampled single cell used in combinations with ^-norm based convex optimization. 

Finally, we applied the developed estimation method to artificial data of a 
medium size bistable system modeling the caspase activation cascade. It could be 
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shown that the proposed method yields good estimation results in case of a setup 
which is realistic in terms of noise and amount of available data. 
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